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ABSTRACT 

The astonishing diversity in the observed planetary population requires theoretical efforts and advances in planet formation theories. 
The use of numerical approaches provides a method to tackle the weaknesses of current models and is an important tool to close gaps 
in poorly constrained areas such as the rapid formation of giant planets in highly evolved systems. So far, most numerical approaches 
make use of Lagrangian-based smoothed-particle hydrodynamics techniques or grid-based 2D axisymmetric simulations. 

We present a new global disk setup to model the first stages of giant planet formation via gravitational instabilities (GI) in 3D with 
the block-structured adaptive mesh refinement (AMR) hydrodynamics code enzo. With this setup, we explore the potential impact of 
AMR techniques on the fragmentation and clumping due to large-scale instabilities using different AMR configurations. Additionally, 
we seek to derive general resolution criteria for global simulations of self-gravitating disks of variable extent. 

We run a grid of simulations with varying AMR settings, including runs with a static grid for comparison. Additionally, we study the 
effects of varying the disk radius. Adopting a simple thermodynamical profile, corresponding to a marginally stable disk ( Q init = 1), 
we validate the numerical robustness of our model for different spatial extensions, from compact to larger, extended disks. The 
physical settings involve disks with R disk = 10,100 and 300 AU, with a mass of M disk « 0.05 M e and a central object of subsolar 
mass (M* = 0.646Af 0 ). To validate our thermodynamical approach we include a set of simulations with a dynamically stable profile 
(Qinu = 3) and similar grid parameters. 

The development of fragmentation and the buildup of distinct clumps in the disk is strongly dependent on the chosen AMR grid 
settings. By combining our findings from the resolution and parameter studies we find a general lower limit criterion to be able to 
resolve GI induced fragmentation features and distinct clumps, which induce turbulence in the disk and seed giant planet formation. 
Irrespective of the physical extension of the disk, topologically disconnected clump features are only resolved if the fragmentation- 
active zone of the disk is resolved with at least 100 cells. The latter corresponds to a minimum requirement for all global disk setups. 
Our simulations illustrate the capabilities of AMR-based modeling techniques for planet formation simulations and underline the 
importance of balanced refinement settings to reproduce fragmenting structures. The clumps in our models are migrating inward and 
are eventually destroyed because of tidal disruptions, reflecting the absence of radiative feedback from the central star, which may 
stabilize the clumps on larger scales. We expect that the inclusion of additional physics, like a radiation transport mechanism and 
the formation of sink particles, will provide a detailed framework to study the formation of planets via gravitational instabilities in a 
global disk view. 
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1. Introduction 


Since the dawn of observational exoplanet detection in re¬ 
cent decades the ever-growing number and variety of planetary 
species gives rise to renewed theoretical efforts for understand¬ 
ing the formation of these bodies (e.g., |Boss|1997[ [Bodenheimer 
|et al.||2000| ). The ultimate goal of achiev ing a consistent ex pla- 
nation for the planetary population (e.g., |Benz et al.||2Q14| ) can 
only be adressed by pushing observational and theoretical tools 
to their limits. 


The two major theories of planet formation, core accretion 
(CA) and formation via gravitational instabilities (GI), have the 
potential to explain a wide variety of planetary species. Whereas 
it is widely accepted that terrestrial planets are formed via a 
bottom-up model like CA (e.g.,|Raymond et al.|2006{|Morbidelli| 


et al.|201 2 ; Raymond et al. 2014), the formation of more massive 


obje cts, like gas giants, can in principle be unders tood in terms of 
CA (Hel led et al. 2014| ) or GI ( [Mayer et al.|2003 1. Recent efforts 
(e.g.^Boley 2009) are dedicated to a synthesis of both models in 
a more general framework of planet formation. Scenarios based 
on GI, however, can be particularly relevant for the formation 
of planets around metal-poor stars, where dust grain coagulation 
becomes increasingly difficult (e.g., Setiawa n’et al.|20lO||John- 


|son et al.]|2012| ), and the buildup of supermassive and Popula- 


tionlll stars (e.g., Inayoshi & Haiman 20 l4[|Latif & Sc hleicher 


|20l5]|Schleicher et al.|2015]) 


In this work we aim for an investigation of the capabilities of 
an alternative technical treatment of GI to achieve a better under¬ 
standing of the buildup of massive planets. This includes gaseous 
planets like Jupiter and Saturn in our own solar system, as well as 
similar or even more massive planets in extrasolar systems (e.g., 
[Borucki et al.|2010| or planets around exotic configurations like 
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in the pre-cataclysmic binary system NN Serpentis. The latter 
system is an example for the class of evolved binaries, for which 


et al. 2013). NN Serpentis recently underwent detailed investi- 

gations (e.g., Haefner et al. 2004; Hessman et al. 2011, 

Horner 

et al. 

2012; Marsh et al. 2014; Volschow et al. 2014; 

Parsons 

et al. 

12014), indicating a second-generation scenario where the 


planets were formed after a common envelope (CE) event. Sub¬ 
sequently, [Schleicher & Dreizler| |2014| ) argued for these plan¬ 
ets to be consistent with the formation via GI and introduced a 
semianalytical model, which further motivates studies in com¬ 
pact self-gravitating disks. In the NN Serpentis scenario, the age 
of the white dwarf of <5 10 6 yr provides an upper limit on the 
timescale for planet formation, therefore favoring GI with re¬ 
spect to CA and indicating a formation scenario in a very com¬ 
pact environment (< 10 AU). 

So far, the hydrodynamical treatment of GI in protoplane¬ 
tary disks mostly relies on a description via smoothed-particle 
hydrody namics (SPH) (e.g.,|Mayer et al.||2007]) or axi symmet- 
ric (e.g., |Boss|2003||Boley et al.|2007 JMeru & Bate|2011| ) static 
grid approaches without adaptive mesh refinement (AMR). Only 
a minority of grid approaches employ a mechanism for adaptive 
and dynamic refinement of the initial mesh (e.g., Paardekooper 
& Mellema 2006 ). In the case of Cartesian grids, the only at¬ 
tempt we know about is the project of |Mayer & Gawryszczak| 
(2008]) and |Gawryszczak et al.| ( |201Q] ), who compared flash and 
gasoline simulations of a self-gravitating disk. In general, the 
use of more than one numerical approach yields the advantage 
of ruling out systematic errors, introduced by the chosen numeri¬ 
cal treatment (e.g.,[Agertz et al. |2007| . Advantages of grid codes 
employing finite volume methods in general are their built-in op¬ 
tion to conserve mass and momentum. Regarding the treatment 
of astrophysical phenomena they offer a precise modeling of tur¬ 
bulence in hydrodynamical instabilities ( [Agertz et al.|2007[ ) and 
a more robust treatment of shocks in comparison with SPH ap¬ 
proaches ( [Price & Federrath|2010| ). Additionally, classical SPH 
appro aches tend to overrate viscosity parameters (Price 2008, 
|2012[ ), which is not the case for grid-based implementations. 
However, recent developments in the SPH community extenu- 
ate many of th ese problems (see, e.g., Read & Hay field 2 012[ 
Hopkins 2013| ). In all numerical codes the chosen geometry dic¬ 
tates the implementation and discretization of the physical equa¬ 
tions. In the case of disk simulations, one could argue that the 
natural geometry of choice would be axisymmetric. However, 
we see some principle advantages in the use of Cartesian ge¬ 
ometry, which motivates comparative studies. For example such 
codes do not include a built-in singularity (i.e., the center of ori¬ 
gin) and handle all computed regions equally. The latter means 
there is no built-in refinement in inner regions, which might 
avoid symmetry-dependent solutions. Including the use of AMR 
schemes enables us to specificially choose regions of greater in¬ 
terest and neglect others. This can be particularly relevant when 
fragmentation occurs because of turbulent motion in the outer 
disk parts and when the deviations from axisymmetry are more 
prominent. 

We use the Eulerian block-structured AMR hybrid code 
enzo (O’She a et al.|2004] Bryan et al.|2014 enzo-project.org), 
which enables massively parallel simulations and supports a 
wide variety of astrophysical problems, i.e., hydrodynamics, 
ideal and non-ideal magnetohydrodynamics, and N-body dy¬ 
namics, including self-gravity for fluids and gas chemistry. The 
fluid flow in enzo is evolved using a finite-volume discretization, 
the ideal gas dynamics calculations are treated with the piece- 
wise parabolic method (PPM), which is a higher-order Godunov 
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scheme ( [Bryan et al.|1995| ) with an accurate piecewise parabolic 
interpolation and a non-linear Riemann solver for shock condi¬ 
tions. The solver is third-order accurate in space and second or¬ 
der accurate in time for fixed time stepping and formally second- 
order accurate for variable time stepping, enzo has been tested on 
a variety of typical fluid flow test problems ( [Bryan et al.|2014| ). 

The code employs the implementation of a (structured) adap¬ 
tive mesh refinement (AMR) scheme, first developed by Berger 


& Colella| ( |1989| ), which enables higher resolution than afford¬ 


able via uniform grids by introducing additional finer mesh 
structures on a coarse uniform grid when necessary (which can 
be defined with a variety of conditions). This scheme utilizes 
an adaptive hierarchy of rectangular grid patches, which cover a 
snippet of space with a certain resolution, refined from the top 
to the bottom of the hierarchy. When the solution evolves and 
interesting regions develop, finer meshes are placed within the 
coarse grid, enabling a higher resolution of these structures and 
not draining too much computing power for less interesting re¬ 
gions. This facilitates a variety of options for modeling GI in 
protoplanetary disks, in which we need to deal with dynamical 
effects on a large range of spatial scales. 

We perform 3D simulations of gravitational instabilities in 
protoplanetary disks with varying disk radius and resolution set¬ 
tings using adaptive mesh refinement techniques. In Section[2]we 
describe the physical and numerical setup of the simulations. In 
Section [3] we compare the initial and evolved state of stable and 
unstable disks and investigate the effect of different refinement 
settings. In Section [4] we verify the setup for various disk radii 
and demonstrate its ability to model systems of a wide range of 
disk extensions. Finally, we summarize our findings in Section 

a 


2. Numerical experiments 

2.1. Numerical setup 

For_a complete description of enzo’s code structure we refer to 
|Bryan etltL] ( |2014| . To provide a sufficient explanation of our 
numerical methods, however, we briefly outline the major nu¬ 
merical recipes used in the calculations. 


Hydrodynamic equations 

In our simulations, we use the piecewise parabolic method 
(PPM), which was originally developed by Colella & Woodward 


(1984) and was adapted for the enzo code by Bry an et al.| (|1995 ) 
in a direct Eulerian fashion. In this implementation, the govern¬ 
ing equations are dimensionally split and rewritten (as example 
in one dimension) in conservative form as 


dp dpv 
dt dx 


dpv dpv 2 dp 
dt + dx + dx 
dpE dpvE 
dt dx 


= P9> 


= pvg , 


( 1 ) 

( 2 ) 

( 3 ) 


in which x and v refer to the one-dimensional position and ve¬ 
locity of the gas, g to the acceleration at the cell center, p to 
the gas density, and E to the total fluid energy density. The one¬ 
dimensional solution to these equations can be found by com¬ 
puting the effective left and right states at each grid boundary. 
This is done by constructing a piecewise parabolic solution of 
density, velocity and total energy and then averaging over the 
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distance that each wave can travel in the specific time step. Then 
a Riemann problem with these effective left and right states is 
solved and all quantities are updated. The standard Riemann 
solver for the PPM method in enzo approximates the appearing 
rarefaction waves, which are traveling left and right, as shocks. 
The solution is then found using an iterative approach. For the 
case of occurring numerical problems due to the Riemann solver 
we used, enzo employs a fallback mechanism. This allows the 
code to switch the Riemann solver to the more diffusive HLL 
solver at particular inte rfaces, where negative de nsities or en¬ 
ergies occur (see, e.g., Lemas ter & Sto ne 2009). Finally, the 
fluxes in between the grid patches are constructed. The full, 
one-dimensional mathematical description and the interpolation 
methods are given in fColella & Woodward] (|1984| and |Bryan| 
[etaLl ( fT995| . 


Because of approximations in this process, |Bryan et ah] ( |2014| ) 
estimates the resolution of the gravitational force to be approxi¬ 
mately twice as coarse as the corresponding refinement level. 


Equation of state 

The total fluid energy density is given by 


E = 



(7) 


with the thermal energy density e governed by the equation of 
state 



( 8 ) 


Dual energy formalism 


As long as the structures in the simulations are well resolved and 
the Mach numbers are on a reasonable scale (i.e., Ma < 100), the 
aforementioned solver works well. However, if hypersonic bulk 
flows with i^therm/^tot ~ 10 -3 appear in the simulation, the nu¬ 
merical situation becomes disastrous: the ratio of kinetic energy 
Ek to gas internal energy e can approach numbers that are much 
too high. Hence, the pressure (proportional to E- E kin ) is the dif¬ 
ference between two extremely large numbers, which would be 
a rather disadvantageous numerical situation with major sources 
for errors, especially for the temperature distribution. To over¬ 
come this problem, enzo separately solves the internal energy 
equation 

de p 

— + v • Vc = — —- V • v, ( 4 ) 

ot p 

with thermal energy density e and pressure p , which is then used 
in hypersonic flows for calculating the pressure and thus the tem¬ 
perature distribution only. The gas dynamics (i.e., the hydrody¬ 
namic routines) should be unaffected to avoid introducing new 
errors. This is achieved by choosing a selection criterion for the 
pressure and thus separating both formalisms via 

I p(y- 1 )(£- v 2 / 2 ), (E-\ 2 /2)/E > T) 1 
p \ p(y- l)e, (£ - v 2 /2)/£ < 77 /• | ; 


If the barrier parameter q is chosen to be small enough, the dy¬ 
namical effect of the dual energy formalism is approximately 


Bryan et al. ( 2014) list 77 = 10 3 as the used standard value, 


zero. 

consistent with usual truncation errors in the simulations. 


Gravity 

Of particular importance for our planet formation model is the 
inclusion of the self-gravity of the gas. To achieve this, Poisson’s 
equation 

V 2 0 = 4/rGp, ( 6 ) 


with the gravitational potential 0 , the gravity constant G, and 
the gas density p, is solved in a multistep process for all cells. 
The cloud-in-cell (CIC) interpolation method from Hockney & 
East wood] ( |1988| ) is used to approximate the gas distribution as 
a spatially-discretized density field, from which the gravitational 
potential is generated via a fast Fourier transform. To get to know 
the acceleration at each (sub)grid a finite difference scheme with 
Dirichlet boundary conditions is used and protected against os¬ 
cillations through buffer zones surrounding the active grid zones. 


This is the equation of state for an ideal gas with the ratio of 
specific heats y. 

Computational domain 

To prevent interactions of the disk material with the boundaries 
of the computational domain, we choose the physical size of the 
box to be (Rdisk • 10) 3 . Thus, there is enough space for the disk to 
evolve and to be embedded in a diffuse medium. The boundaries 
of the domain behave as reflecting walls, i.e., 

q(-x) = q(+x), (9) 

with an arbitrary field quantity q , and with the velocity perpen¬ 
dicular to the boundary direction reversed 


v x (-x) = -v x (+x). 


( 10 ) 


Refinement criteria 


The refinement to the next (deeper) level of the grid patches in 
our simulations depends on the corresponding cell mass and a 
so-called Jeans length criterion. The cell mass refinement works 
such that it mimics a Lagrangian method in trying to keep a fixed 
mass resolution. Thus, if the mass of a specific cell is 

Mg = p(Ax) d > pflag(Ax root ) d r e ‘ l = Mjupiter/ 1000 , ( 11 ) 


the cell is refined to a deeper level (if allowed by the maximum 
refinement level). Here, pfi ag is the required density on the root 
grid, Ax root is the root grid cell spacing, r the refinement fac¬ 
tor (usually 2 ), l the level and 6 / is an aggressiveness parame¬ 
ter, to make the refinement super-Lagrangian ( 6 / < 0) or sub- 
Lagrangian ( 6 / > 0). 

The second approach refines the Jeans length 


4 / 



( 12 ) 


by a fixed number of cells (based on Truelove et al. [1998| ) if 


( t±bt y /2 

{N}Gpm H ) 


(13) 


with Nj = 32 the required number of cells per Jeans length. 
This is especially valuable since we deal with fragmentation ef¬ 
fects, where contraction plays a key role in the later building of 
clumps. 
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2.2. Initial conditions 


3D density profile 

To begin with, we set the column density structure of the disk 
according to a classical Mestel power-law profile ( |Mestel|1963| ) 

n out 


2(r) = Xo- 


(14) 


with r out denoting the outer radius of the disk and £ 0 th e column 
density at this radius. Since the contribution to mass and angular 
momentum at small disk radii can be neglected we assume that 
the disk extends to r = 0. This leads to a normalization of the 
column density as 

1 A/disk 


2n 


(15) 


with the overall disk mass chosen to be M^k ~ 0.05 M©. Using 
hydrostatic balance the analytical solution to resolve the Gaus¬ 
sian density distribution in the vertical direction is 


p(r, z) ■ 


£(r) 


exp 


H(r) V2tt \ 2H(r)- 


(16) 


Here, z denotes the distance to the midplane of the disk and H 
the scale height, which is evaluated as ( |Lodato|2008 ) 


*(0 = ^ 
2 \H sg 



■■ 0.627/*, 


(17) 


with the scale height for self-gravitating disks 

7Z sg = djnGY, (18) 

and the scale height for Keplerian disks 

H+ = c s /n. (19) 

In the equation above, the sound speed is given via the Newton- 
Laplace equation 


P 

r • 


( 20 ) 


and th e an gular velocity fl (see next subitem). As mentioned 
by |Lodato| ( f2008| ), this approximation holds for a comparable 
gravitational influence of the disk and the central object, H sg « 
H+. 


Central object 

We choose the host star of the system to be of subsolar mass 

M* = 0.646M 0 , (21) 

which is represented in the simulations as a point mass parti¬ 
cle without spatial extension (denoted in the later figures with 
a single black dot in the middle) and is included in the gravity 
calculations. 


Orbital velocity profile 

Including the disk’s self-gravity, the azimuthal velocity satisfies 
( |Lodato|2008[ ) 

” 2 Id ;■> GM* „ _ Id p 

y y ( 22 ) 


dr 


p dr 


+ 2 7tG2 + 


p dr ’ 


with d> tot the gravitational potential of central object and disk 
material. 
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Thermal profile 


Finally, we need to address the thermal profile of the disk. We 
characterize the disk through the Toomre parameter (Toomre 


C S K c s £l 

7tGS 7tG£ ’ 


(23) 


with Boltzmann’s constant k s . Fixing Q to a constant value then 
translates into a condition for the sound speed and therefore the 
temperature of the gas. The epicyclic frequency k can be approx¬ 
imated with the angular velocity fl = v^/r, when the velocity 
profile throughout the disk is mostly dominated by the Keple¬ 
rian velocity ( |Lodato|2008| . We assume a gas with composition 
similar to the solar system, i.e., the mean molecular weight p is 
set to 2.4 ( [Mayer et al.|2007| ). From this assumption we infer the 
temperature profile using Equation [20] to be 


T(r) = 


yh 


(24) 


The ratio of specific heats is fixed to be y = 1.0001 « 1, accord¬ 
ing to an isothermal equation of state (EOS). 


3. Resolution study 

In this chapter, we investigate the imprint of the refinement level 
on the physical conditions within the disk. We benchmark the 
simulations for disks with Rdisk = 100 AU for the initial condi¬ 
tions and at t - 1.5 Tdisk, i.e., 1.5 orbital times of the most outer 
radius of the disk. An overview of all of our featured simulations 
is given below in Table [I] 


Rdisk [AU] Qinit Pc imax 


Run 

R10i64r2 

R10i64r3 

R10i64r4 

RlQil28r2 

R10sg256 

RlO0i64r2 

RlO0i64r3 

RlO0i64r4 

R100il28r2 

R100sg256 

R300i64r2 

R300i64r3 

R300i64r4 

R300il28r2 

R300sg256 

Qi64r2 

Q3i64r3 

Q3i64r4 

Q3il28r2 

Q3sg256 


10 1 

10 1 

10 1 

10 1 

10 1 

100 1 

100 1 

100 1 

100 1 

100 1 

300 1 

300 1 

300 1 

300 1 

300 1 

100 3 

100 3 

100 3 

100 3 

100 3 


6U 2 

64 3 3 

64 3 4 

128 3 2 

256 3 0 

64 3 2 

64 3 3 

64 3 4 

128 3 2 

256 3 0 

64 3 2 

64 3 3 

64 3 4 

128 3 2 

256 3 0 

64 3 2 

64 3 3 

64 3 4 

128 3 2 

256 3 0 


Table 1: Overview of the settings of all performed simulation 
runs, with the disk radius Rdisk, the Toomre parameter Qi nit , the 
initial number of cells of the coarsest grid in the simulation box 
g c , and the maximum allowed refinement level l max . 
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0 20 40 60 80 100 120 


r [AU] 


Fig. 1: Radial profiles for initial column density X and rotational 
velocity v rot for simulations with Rdisk = 100. The parameter 
g c indicates the resolution of the coarsest grid, l max the deep¬ 
est level of the resolution by dynamic refinement. The continous 
lines show simulations with Q init = 1, the dashed lines simula¬ 
tions with Q init = 3. The density distribution converges toward 
the analytic solution for higher refinement levels, i.e., the density 
gap in the inner part of the disk is smaller for higher refinement 
levels. Additionally, it is smaller for the simulation runs with 
Qinit = 1, which is due to the scale height dependency on the 
Q parameter (compare Equation [23] and Figure [3]). The velocity, 
only nonzero within the disk matter, behaves analogous to the 
density. 


3.1. Initial state 

Figure[l]and[2]show comparisons of the radial profiles of the val¬ 
ues of X, v rot , T and Q of disks with Rdisk = 100 AU and Q m i t = 1, 
Qinit = 3, respectively. All of them converge for higher maxi¬ 
mum refinement levels toward the analytical solution, as given 
in Section [221 Some differences between the simulations with 
Qinit = 1 and Q^n = 3 arise because of the change of scale 
height H by a change in g, according to Equation [15] It is im¬ 
portant to say that the simulations with the same effective resolu¬ 
tion behave very similarly, despite any differences in resolution 
for more coarse grids. Therefore, we restrict the discussion on 
the simulations with g c = 64 3 , while analogous simulations with 
the same effective resolution can be identified for the runs with 
g c - 128 3 . In the following we discuss and explain the differ¬ 
ences in detail. 


Fig. 2: Radial profiles for initial temperature T and effective 
Toomre Q parameter for simulations with Rdisk = 100. The tem¬ 
perature increases outward because of its dependence on Q (see 
Equation [24]). Q is approximated correctly in the inner parts of 
the disk. In outer parts, where the density decreases for numeri¬ 
cal reasons, Q rises, and thus artificially stabilizes the disk mat¬ 
ter. 


The profiles for the surface density X yield two conclusions. 
First, higher refinement means that the disk is much better re¬ 
solved, especially in the inner part. Second, it displays the de¬ 
pendency of H on Q. This is illustrated in Figure [3] which shows 
slices of the initial density profile in a cut through the vertical 
profile of the disk, comparing disks with Rdisk = 100 AU and 
Qinit = 1, Qinit ~ 3, respectively. The different refinement levels 
are indicated with gray boxes. These show that the disk midplane 
is better resolved than the outer parts, thus enabling a higher res¬ 
olution in the regimes important for fragmentation. The vertical 
disk extension is no longer resolved when the scale height drops 
below the cell size at small radii. 


All disks represent the analytical solution for the rotational 
velocity very well, with only minor differences. They are better 
resolved (i.e., better converge toward the analytic solution in the 
inner regions) for a higher refinement level. This is a direct con¬ 
sequence from the density distribution (figures [1] and [3}, since 
the rotational velocity is only defined in the disk materia . 

The temperature profiles show fundamental differences for 
the Qinit = 3 and Qi n n = 1 cases as a result of the Q depen¬ 
dence of the thermal profile (Equation [24]). Fixing the Toomre 
parameter to Q m it = 1 results in lower temperatures throughout 
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Fig. 3: The initial vertical density profile 2, comparing disks with 7?^ = 100 AU and = 3 (top), = 1 (bottom), and 
increasing maximum refinement level from 2-4 (left to right), respectively. The scale height H, according to Equation [T9J changes 
with Q, which is why the Q in i t = 1 disks show a flatter profile. The gray boxes indicate the deepness of the refinement level, from 
level l = 0 outside of the boxes to level / = l max in the most inner box. 


the disk and thus an increased chance for developing fragmen¬ 
tation. As a result, the temperature distributions show a posi¬ 
tive slope with increasing radius, which is adopted to explore 
disks in a marginally stable state. Outside of the disk the temper¬ 
ature rapidly increases to the temperature of the background gas 

T background — 200 K. 

Finally, the effective Toomre Q distributions show the nu¬ 
merical imprint of Q init on the disk as a consequence of dis¬ 
cretization effects. Thus, Q drastically increases in the inner and 
outer parts of the disk, where the density per cell decreases for 
numerical reasons. Despite these limitations Q stays constant in 
the dynamically relevant disk parts, where the density is high 
enough to support fragmentation effects. 


3.2. Evolved state 

A potential weakness of Cartesian grid codes is angular momen¬ 
tum dissipation, which is especially relevant for modeling sys¬ 
tems with high rotational velocities. Therefore, to demonstrate 
the numerical robustness of our setup, we show the evolution of 
the total angular momentum L tota i in the whole simulation box 
in Figure]?] for simulation times up to t = 1.5 Tdisk- 

Overall, the deviations for all simulations go at maximum 
up to +2% of the initial angular momentum. For the simulations 
with Q init = 3, the angular momentum shows a very stable con¬ 
figuration and all simulations show the same general trend. How¬ 
ever, the simulations with Qi n a - 1 differ from each other, de¬ 
pending on the deepest refinement level in the simulation. Here, 
the simulations with the highest resolutions (i.e., with deepest 
refinement level l max = 3 or l max = 4) show a decrease of up to 
-1.5% total angular momentum. 

Now we focus on the imprint of changes in resolution on 
the evolved state of the disks for t - 1.5 Tdisk to spot deviations 
for different refinement depths. Figure [5] visualizes the density 
distribution for a face-on view of the disks for various refine¬ 
ments for Qi n it - 3 and Q in i t = 1. The disks with Q init = 3 do not 
differ qualitatively with increasing refinement level, but show in¬ 
creased resolution and details in the gas streams. 

However, the differences for the Qi n u = 1 disks are dramatic. 
Whereas the disk with l max = 2 does not show any signs of 
clumping or spiral arms, these characteristics start to appear for 
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Fig. 4: Evolution of angular momentum for all simulations with 
Rdisk - 100 and up to t = 1.5 T disk . The maximum deviations 
from the initial state reach up to 2%. 


lmax = 3 and are sharply defined in the simulation with l max = 4. 
The differences due to l max do not only change the overall density 
distribution of the disk. In addition, the planetary (clump) wake, 
associated with each forming clump, differs dramatically from 
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Fig. 5: Face-on column density profiles for simulations with g c = 64 3 and t = 1.5 T^sk- The top row shows simulations with 
Qinit = 3, the bottom row simulations with Q init = 1. The maximum refinement level l max is increasing from left to right from 2-4. 



Fig. 6: Edge-on volume density profiles for simulations with g c = 64 3 and t = 1.5 The top row shows simulations with 
Q init = 3, the bottom row simulations with Q init = 1. The maximum refinement level l max is increasing from left to right from 2-4. 


l max = 3 to l max = 4. For l max - 3, the clump only marginally dis¬ 
turbs the large-scale structure of the spiral arms, forming regions 
of higher density at approximately 50 and 80-90 AU separation. 
In contrast, the simulation with l max = 4 shows very distinct 
spots, where fragmentation occurs. These clumps massively dis¬ 
turb the spiral arm structure, caused by the turbulent motion in 
the fragmenting areas. 


The edge-on views of the density distribution of the disks 
(Figure [6]) reveal structural differences in the vertical gas dis¬ 
tribution with increasing refinement level. The Qi n u = 3 disks 
again show relatively weak differences, however, the inner part 
of the disk flattens with increasing refinement and thus the flared 
disk profile is better resolved. 

The Qi n i t = 1 disks show this flattening as well. Moreover, 
the density peaks in the midplane parts of the disks are more pro- 
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Fig. 7: Radial profiles of column density £ and rotational veloc¬ 
ity v rot for all simulations with Rdisk = 100 and t - 1.5 Tdisk- 
The runs with Q init = 3 approximately follow the initially cho¬ 
sen density distribution. The Q in d = 1 disks show peaks at radii 
with distinct clumps or spiral structure (compare figures [5] and 
[6]). The velocity profiles show only minor deviations from the 
general trend, which are due to the turbluent velocity in the most 
unstable regions. 


nounced with higher refinement level. Comparing the edge-on 
slices of the disks with the face-on projections demonstrates the 
capabilities of three-dimensional simulations. In areas of strong 
fragmentation the disk’s vertical extension is lower compared to 
more stable regions. Thus, the flared profile of the disk is dis¬ 
turbed and shows strong deviations from vertical and axisym- 
metry. 

The visual observations are supported by the quantitative 
analysis of the radial profiles in figures [7] and [5] 

Figure [7] shows the density distributions of the evolved state. 
The simulations with Q init = 3 roughly follow the initally cho¬ 
sen density profile. In the inner and outer parts the profile is 
marginally spread out. The increase of density in the inner parts 
is the result of our analytic velocity approximation. Therein we 
assume the disk to be extended all the way down to the center. 
However, since the disk’s inner extension is cut off at some ra¬ 
dius (see Section HU the pressure gradient at the inner radius 
of the density distribution is too high. Therefore, the effective 
velocities in the inner part of the numerical implementation are 
too small. Thus, some of the gas on lower orbits rapidly migrates 
toward the center until an equilibrium state is reached. A similar 


Fig. 8: Radial profiles of the temperature distribution and effec¬ 
tive Toomre Q parameter for all simulations with Rdisk = 100 and 
t = 1.5 Tdisk’ As the density distribution the temperature shows a 
spread-out for the inner and outer parts. This is much more pro¬ 
nounced for the simulations with Q init = 1, whereas the runs with 
Qinit = 3 roughly follow the initial distribution. For Q init = 3, the 
effective Q parameter stays well above the threshold for marginal 
stability. In comparison, the effective Q in the Qimr = 1 runs dips 
below the threshold at the positions with increased density (i.e., 
fragmenting regions, compare Figure [7]). 


effect occurs in the outer parts, where the gas is spread out until 
a smooth transition from high to low density is reached. 

In comparison, the density distributions of disks with Q init - 
1 show deviations from the initially chosen profile. The turbulent 
and self-gravitating disk material induces stochastic fluctuations 
in the velocity. This, in turn, initiates the development of spiral 
arms and fragmentation, which is reflected in the radial profiles 
as peaks. Overall, the most massive clump can be seen for the 
simulation with l max = 4 at ~ 52 AU, reaching column densities 
up to 10 3 -^2 • 

r cm z 

As discussed above, the velocity profiles for all disks reach 
an equilibrium state, which follows the analytical solution in 
most parts of the disk. The small deviations in the distribution for 
the 2-1 mns emerge from the spiral arms and fragmentation, 
which are more pronounced for higher maximum refinement. 

The temperature distribution in Figure [8] behaves analo¬ 
gously to the density in Figure [7] Thus, it features a spread in 
the inner and outer parts, more remarkable for lower resolutions. 
Higher resolutions however are able to sustain the initial profile 
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Fig. 9: Number of lowest-level clumps for all simulations and 
t - 1.5 Tdi S k. The trend goes for just a few clumps for low maxi¬ 
mum resolutions to many clumps for high maximum resolution. 
Different combinations of g c and l max but comparable resolution 
on the lowest level (i.e., green/cyan and blue/magenta) show ap¬ 
proximately equal numbers of lowest-level clumps. 


better, which is true for both Q init = 3 and Qmif = 1 simulations. 
All of them reflect the spread of the gas and feature an approxi¬ 
mately similar temperature distribution to the outer edge. 

The inital Q state is well preserved for the gmzY = 3 runs. 
However, the strong instabilities and the resulting fragmentation 
in the simulations with Q init = 1 lead to minima in the effective 
Q below the unity threshold where the density reaches its highest 
values. 

As an attempt to quantify the behavior of the fragmenta¬ 
tion for different refinement settings, we use the clump finding 
method of |Smith et~aL] ( |2010| to detect topologically discon¬ 
nected structures within the dataset. Its principle mechanism is to 
create a single contour in between p min and p max over the whole 
computational domain. The lower value is then continously in¬ 
creased, until it reaches the maximum value. During that process 
disconnected structures are identified as separate contours and 
are treated as individual structures in which the routine contin¬ 
ues recursively. To get an idea of the smallest structures in our 
simulations, we only print the far-end leaves (“clumps”) of this 
treelike structure, i.e., topologically disconnected regions featur¬ 
ing approximately the highest density in their specific contours. 

Figure [9] illustrates our findings, showing a histogram with 
statistics for all simulations with Q init = 1. The number of 
clumps increases with the maximum level of refinement, i.e., we 
find more clumps for higher resolution. This underlines the im¬ 
portance of a thorough resolution to achieve correct structures 
throughout the grid. 


4. Parameter study 

In this section we study the influence of the disk extension on 
the resolution dependent outcome of our model. In general, ob¬ 
servations of protoplanetary disks show a wide variety of disk 
morphologies (e.g., |Avenhaus et al.|[20l4{ |Garufi et~a l. 2014). 
Therefore, disk characteristics are likely to vary in radial extent. 
Additionally, some exotic systems like NN Ser are expected to 
host very compact systems of only a few AU radius ( [Schleicher] 




Fig. 10: Radial profiles at t = 1.5 Tdisk for column density £ 
and rotational velocity v rot for all simulations with the highest 
resolution on the finest grid, normalized in radial direction. The 
simulation with Q init = 3 is added for comparison. The density 
distributions for all radii show similarities. However, the simu¬ 
lation with r = 100 AU and Q init = 1 shows its single peak 
at ~ 52 AU and the negative slope at the outer edge of the disk 
is much shallower than the slope of the other simulations with 
Q init = 1. The velocity distributions for all radii show a similar 
trend. The simulation with r - 10 AU abruptly goes down in 
velocity. This behavior might be due to the very sharp transition 
from disk material (with approximately Keplerian velocity pro¬ 
file) to surrounding material (with approximately zero velocity). 


& Dreizler|2014| . Thus, we test our setup on its ability to model 
disks of varying radius. Note that all our simulated disks have 
the same overall mass Mdisk and because of the Q dependency 
of our model, the temperature is adjusted to be able to reproduce 
global Q unstable configurations. 

Figure [TO] and [TT] show radial profiles for the parameters £, 
v rot , T and effective Q for all simulations with the highest resolu¬ 
tion, i.e., with maximum refinement level l max = 4. In principle, 
the simulations with different initial disk radii reproduce similar 
results and the distributions of the various physical parameters 
show similar trends. However, there are some differences, which 
mainly arise because of the single peak in the density distribution 
of the simulation with Rdisk = 100 AU and Qmfr = 1 at around 
52 AU separation from the central object. This peak, resulting 
from a vary massive clump (compare Figure [l2|), influences the 
whole disk structure, and drains material from the other parts 
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Fig. 11: Radial profiles at t = 1.5 Tdisk for temperature T and 
effective Toomre Q parameter for all simulations with the highest 
resolution on the finest grid, normalized in radial direction. The 
simulation with Q init = 3 is added for comparison. Again, we see 
a different behavior of the simulation with Rdisk = 100 AU and 
Q init = 1, displaying a smaller increase in temperature because 
of its smoother transition zone. The dips where the Q parameter 
goes below 1 are consistent with the density peaks in figures [TO] 
and 


12 


of the disk, which results in a different average distribution of 
matter in the disk, compared to the simulations with Rdisk = 10 
AU and Rdisk = 300 AU. While we expect this peak essentially 
to represent the stochastic nature of fragmentation, its presence 
leads to differences for many other parameters as well. For exam¬ 
ple, the velocity decrease is shallower, especially in comparison 
with the Rdisk = 10 AU run, which shows a sharp transition at 
its outer edge from disk material to surrounding material. Ad- 
dtionally, the temperature distribution shows an approximately 
constant distribution up to 1.3 Rdisk, because the (cold) disk ma¬ 
terial (in comparison with the surrounding gas) is extended over 
a larger radius and therefore the transition from disk gas to sur¬ 
rounding gas happens at larger radii. 

The Toomre Q radial profiles can be matched very well with 
the profiles for density and temperature. Thus, the most pro¬ 
nounced dip in the Q radial profile is found at approximately 
0.5 Rdisk separation for the simulation with Rdisk = 100 AU. The 
simulations with Rdisk = 10 AU and Rdisk = 300 AU show a very 
similar (radially) normalized behavior. Their most pronounced 
Q dips are at separations 0.5 Rdisk and 0.8 - 0.9 Rdisk- Matching 
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this observation with Figure [12] we see that the highest amount 
of fragmentation can be spotted in these regions. 

Figure [l2| shows a final impression of the evolved structures 
for all simulations. The turbulent structures of the disk can be 
seen in the column density, as well as in the volume density plots. 
The projected face-on views provide detailed insight into the 
fragmented disk structures. The edge-on views represent very 
well the flared disk structures and reveal the regions of the high¬ 
est density, where the most massive clumps are located. 


5. Conclusions 


We performed simulations of self-gravitating, massive proto¬ 
planetary disks using the AMR code enzo. Our physical setup 
was m otivated by the semianalytical approach of Schleicher & 
|Dreizler| ( |2014| ) for modeling the characteristics of a formation 
via GI in diverse planet formation environments. We calculated 
the time evolution of disk configurations close to the Qi n u = 1 
threshold with disk radii of Rdisk = 10,100 and 300 AU, varying 
grid settings, and the time evolution of a stable disk (Qinit = 3) 
with Rdisk =100 AU to validate our approach. 

The Qi n it = 1 disks display the onset of large-scale GI on the 
orbital timescale, yielding gas fragmentation and the formation 
of clumps throughout the disk material. The imprint of AMR 
effects plays a major role in correctly modeling the fragmen¬ 
tation in the disk. Only the highest resolution levels yield pro¬ 
nounced signals of fragmentation and single peaks in the devel¬ 
oping spiral arms. Additionally, the structure of the spiral arms 
and the fragmentation process is qualitatively different for dif¬ 
ferent refinement levels. Whereas for low resolution the spiral 
arms are less pronounced and the forming clumps differ only 
slightly from their surrounding medium, with higher resolution 
the clumps build sharp column density peaks. 

To summarize these findings, we are able to define a lower 
limit for the resolution in global AMR simulations of GI induced 
fragmentation in self-gravitating protoplanetary disks. A com¬ 
parison of the imprint of the maximum refinement level on the 
fragmentation structure and the building of clumps (figures [5]and 
[9]) yields a minimum ratio of fragmenting disk radius Rdisk to 
resolution level r/ y / (in physical units). Thus, to induce the devel¬ 
opment of spiral arms in the disk this ratio has to be 


Rdisk/rivi ^ 50, 


(25) 


corresponding to our simulations with l max = 3. However, to 
be able to reveal further fragmentation and the building of dis¬ 
tinct clumps in the disk spiral structures and resolve the planetary 
wake around them, one has to satisfy the more rigid criterion of 


Rdisk/rivl ^ 100, 


(26) 


which corresponds to our simulations with l max = 4. These crit- 
era, giving the number of horizontal cells the disk radius should 
contain, should be understood as a lower limit to resolve frag¬ 
mentation structures in a global disk view. Running higher res¬ 
olutions will most likely yield even better and more resolved 
structures. 

As another valuable feature of the 3D implementation of 
our setup, disks of all radii conserve their flared structure very 
well and resolve deviations in scale height during the fragmen¬ 
tation process. This might be especially important to investigate 
deviations from vertical and axisymmetry in the planet forma¬ 
tion process. This is in agreement with the findings of Ma yer &| 
Gawryszczak (2008), who underline the importance of (initial) 
resolution in the vertical direction to resolve clumping effects. 
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Fig. 12: Top row : Density projections of all disks with Q init = 1 and l max = 4, t = 1.5 T disk . Bottom row : Vertical density slices of all 
disks with Q init = 1, l max = 4 and t = 1.5 T disk . 


Even if it is in general hard to compare the outcome of SPH cal¬ 
culations with those from grid-based models, our resolution de¬ 
pendent results are similar to those of |Nelson1 ( [2006j ) concerning 
the importance of vertical resolution in SPH simulations. 


The total angular momentum in all simulations is approxi¬ 
mately conserved (with total errors of up to 2%) during the time 
evolution of the simulations. The marginal loss in angular mo¬ 
mentum, shown in Figure [4] can be explained with discretiza¬ 
tion effects due to the grid structure of the simulations. From a 
numerical point of view the simulations presented here can be 
understood as a first step toward a more realistic coverage of 
turbulence effects in accretion disk scenarios with high Mach 
numbers. As found by Federrath et al. ( [201 lj ) the energy injec¬ 
tion scale of gravity-driven turbulence is close to the local Jeans 
scale. Therefore, in AMR simulations in which gravitational en¬ 
ergy is converted into turbulent energy, it is crucial to resolve the 
Jeans length by (at least) 32 ( [Federrath et al.||20TTj ) or even 64 
( [Latif et al.|2013a|b| ) cells. Therefore, the maximum refinement 
level chosen in our simulations is a constraint for even better 
angular momentum conservation. 


Additional ideas to speed up simulations of self-gravitating 
accretion disks and to enable simulations on a wider range of 
spatial scales are necessary to eventually cover the range from 
the global disk view down to smaller scales as, e.g., the physics 
of circumplanetary disks (Morbidel li et al.|[2013| ). A first step 
toward this goal might be the introduction of a more specific 
refinement criterion, which only covers areas that can be asso¬ 
ciated with long-lasting clumps, which might eventually end up 
on a planetary mass scale. 


A particular interesting feature of the simulations from a 
physical point of view is the formation and rapid inward mi¬ 
gration of clumps. They migrate inward within a few orbital 
times, which is roughly consistent with the timescale of Type 
I migration ( [Helled et al.|2014| ). Apart from the relatively small 
masses, the rapid inward migration is consistent with the find¬ 
ings of |Michael et al.| ( [2011| ) and |Baruteau et al.| ( |2011| ). The 
migration in unstable and turbulent disks is driven by several 
forces, such as clump-disk interactions (like Lindblad torques) 
or stochastic torques. Identifying the major driving forces in our 
scenario and to find a stabilization mechanism for the migrating 
clumps are subjects for further research. 


To overcome these fast migration scenarios, [Meru et aL| 
( 2014[ ) proposed the sculpting of long-lasting dust rings through 
pressure maxima in the disk. However, a scenario to explain the 
migration stalling purely based on the thermal evolution of the 
gas might be the influence of radiative feedback from the cen¬ 
tral star. Following Chiang & Goldreich ( |1997| ), the photoheat¬ 
ing would transform a significant p art of the disk to be stabilized 
against fragmentation instabilities ( [Schleicher & Dreizler|2014| ). 
Clumps migrating toward such stabilized regions might still be 
able to grow and evacuate the disk in this region. This stalling 
of the inward migration and possible opening of gaps in the disk 
would prevent them from being destroyed by tidal interactions 
( |Zhu et al.|20l2] ). 


The timescale of the disk’s evolution, comparable to the 
orbital timescale, and the fast fragmentation and formation of 
clumps is of major interest for the formation of planets in exotic 
systems like the NN Serpentis binary system. Since its evolution¬ 
ary timescale is determined by the cooling age of its white dwarf 
(Hessman et aL| |2011| ), the formation of the planets must hap¬ 
pen on a timescale of < 10 6 yr. Thus, the rapid fragmentation 
and clumping opens up a possible way to create the proposed 
NN Ser planets in the appropriate time. However, further obser¬ 
vational coverage and more realistic theoretical investigations of 
this system are needed. 

We have presented here a systematic study exploring the GI 
in self-gravitating protoplanetary disks with AMR simulations. 
We expect that this technique can provide additional insight 
into the formation of massive self-gravitating clumps in future 
simulations, which will complement the existing numerical ap¬ 
proaches well. This technique can be combined with additional 
physics modules like the chemistry pa ckage krome ([Grassi et al. 
2014), radiation transp ort techniques ([Wise & Abel|2011| ) and a 
sink particle algorithm ( [Latif et al.|2013c| ) for an improved mod¬ 
eling of the formation of planets. The simulations presented here 
will be particularly valuable as a reference model to understand 
the evolution of the GI in the absence of additional processes 
that can complicate the situation, but already shed light on the 
important role and necessity of high resolution in grid-based cal¬ 
culations. 
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